// C++ source file Sash.cpp
// Implementation of the SASH index for approximate similarity search,
// as described in 
//   Michael E. Houle (author),
//   "SASH: a Spatial Approximation Sample Hierarchy for Similarity Search",
//   IBM Tokyo Research Laboratory Technical Report RT-0517, 5 March 2003.
// and
//   Michael E. Houle and Jun Sakuma (authors),
//   "Fast Approximate Search in Extremely High-Dimensional Data Sets",
//   in Proc. 21st International Conference on Data Engineering (ICDE 2005),
//   Tokyo, Japan, April 2005, pp. 619-630.
//
// Copyright (C) 2004-2006, Michael E. Houle,
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//
// 1. The source code and derived binary forms may be used only for
//    non-commercial, non-profit research purposes.
//
// 2. Redistributions of source code must retain the above copyright
//    notice, these conditions, and the following disclaimer.
//
// 3. Redistributions in binary form must reproduce the above copyright
//    notice, these conditions, and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 4. The names of its contributors may not be used to endorse or promote
//    products derived from this software without specific prior written
//    permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
// OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Comments, bug fixes, etc welcome!
// Contact e-mail address: meh@nii.ac.jp, meh@acm.org


/**
 * Sash.cpp:    SASH index for approximate similarity search.
 *
 * Author:      Michael E. Houle
 * Date:        4 Jan 2006
 * Version:     1.0
 */

#include "Sash.h"

////////////////////////////////////////////////////////////////////////
//                                Sash                                //
////////////////////////////////////////////////////////////////////////


  //////////////////////////////////////////////////////////////////////
  //                         Public Methods                           //
  //////////////////////////////////////////////////////////////////////


  /**
   * Constructor using default seed value for
   * random number generator initialization.
   */

  Sash:: Sash ()
  //
  {
    data = NULL;
    size = 0;
    maxParents = 4;
    maxChildren = 16;
    internToExternMapping = NULL;
    sampleSizeList = NULL;
    levels = 0;
    parentIndexLList = NULL;
    parentDistLList = NULL;
    parentLSizeList = NULL;
    childIndexLList = NULL;
    childDistLList = NULL;
    childLSizeList = NULL;
    query = NULL;
    distFromQueryList = NULL;
    storedDistIndexList = NULL;
    numStoredDists = 0;
    numDistComps = 0UL;
    levelQuotaList = NULL;
    queryResultIndexList = NULL;
    queryResultDistList = NULL;
    queryResultSize = 0;
    queryResultSampleSize = 0;
    scratchIndexList = NULL;
    scratchDistList = NULL;
    scratchListSize = 0;
    verbosity = 0;
    numOrphans = 0;

    stringBuf = new char [SASH_BUFSIZE_];
    stringBuf[0] = '\0';

    genInt.seed (314159UL);
    seed = 314159UL;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Constructor using seed for random number generator initialization.
   */


  Sash:: Sash (unsigned long seed)
  //
  {
    data = NULL;
    size = 0;
    maxParents = 4;
    maxChildren = 16;
    internToExternMapping = NULL;
    sampleSizeList = NULL;
    levels = 0;
    parentIndexLList = NULL;
    parentDistLList = NULL;
    parentLSizeList = NULL;
    childIndexLList = NULL;
    childDistLList = NULL;
    childLSizeList = NULL;
    query = NULL;
    distFromQueryList = NULL;
    storedDistIndexList = NULL;
    numStoredDists = 0;
    numDistComps = 0UL;
    levelQuotaList = NULL;
    queryResultIndexList = NULL;
    queryResultDistList = NULL;
    queryResultSize = 0;
    queryResultSampleSize = 0;
    scratchIndexList = NULL;
    scratchDistList = NULL;
    scratchListSize = 0;
    verbosity = 0;
    numOrphans = 0;

    stringBuf = new char [SASH_BUFSIZE_];
    stringBuf[0] = '\0';

    genInt.seed (seed);
    this->seed = seed;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Destructor.
   */

  Sash:: ~Sash ()
  //
  {
    int i;

    data = NULL;

    if (internToExternMapping != NULL)
    {
      delete [] internToExternMapping;
      internToExternMapping = NULL;
    }

    if (sampleSizeList != NULL)
    {
      delete [] sampleSizeList;
      sampleSizeList = NULL;
    }

    if (parentIndexLList != NULL)
    {
      for (i=0; i<size; i++)
      {
        if (parentIndexLList[i] != NULL)
        {
          delete [] parentIndexLList[i];
          parentIndexLList[i] = NULL;
        }
      }

      delete [] parentIndexLList;
      parentIndexLList = NULL;
    }

    if (parentDistLList != NULL)
    {
      for (i=0; i<size; i++)
      {
        if (parentDistLList[i] != NULL)
        {
          delete [] parentDistLList[i];
          parentDistLList[i] = NULL;
        }
      }

      delete [] parentDistLList;
      parentDistLList = NULL;
    }

    if (parentLSizeList != NULL)
    {
      delete [] parentLSizeList;
      parentLSizeList = NULL;
    }

    if (childIndexLList != NULL)
    {
      for (i=0; i<size; i++)
      {
        if (childIndexLList[i] != NULL)
        {
          delete [] childIndexLList[i];
          childIndexLList[i] = NULL;
        }
      }

      delete [] childIndexLList;
      childIndexLList = NULL;
    }

    if (childDistLList != NULL)
    {
      for (i=0; i<size; i++)
      {
        if (childDistLList[i] != NULL)
        {
          delete [] childDistLList[i];
          childDistLList[i] = NULL;
        }
      }

      delete [] childDistLList;
      childDistLList = NULL;
    }

    if (childLSizeList != NULL)
    {
      delete [] childLSizeList;
      childLSizeList = NULL;
    }

    query = NULL;

    if (distFromQueryList != NULL)
    {
      delete [] distFromQueryList;
      distFromQueryList = NULL;
    }

    if (storedDistIndexList != NULL)
    {
      delete [] storedDistIndexList;
      storedDistIndexList = NULL;
    }

    if (levelQuotaList != NULL)
    {
      delete [] levelQuotaList;
      levelQuotaList = NULL;
    }

    if (queryResultIndexList != NULL)
    {
      delete [] queryResultIndexList;
      queryResultIndexList = NULL;
    }

    if (queryResultDistList != NULL)
    {
      delete [] queryResultDistList;
      queryResultDistList = NULL;
    }

    if (scratchIndexList != NULL)
    {
      delete [] scratchIndexList;
      scratchIndexList = NULL;
    }

    if (scratchDistList != NULL)
    {
      delete [] scratchDistList;
      scratchDistList = NULL;
    }

    if (stringBuf != NULL)
    {
      delete [] stringBuf;
      stringBuf = NULL;
    }
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Constructs the SASH from an array of data items.
   * A default capacity of 4 parents per node is assumed.
   * The number of items comprising the SASH is returned.
   * The number of pairwise distance computations performed
   *   can be obtained via a call to getResultDistComps.
   */

  int Sash:: build (DistData** inputData, int numItems)
  //
  {
    return build (inputData, numItems, 4);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Constructs the SASH from an array of data items.
   * The maximum number of parents per node must be specified.
   * If a value smaller than 3 is provided, 3 is used instead.
   * The number of items comprising the SASH is returned.
   * The number of pairwise distance computations performed
   *   can be obtained via a call to getResultDistComps.
   */

  int Sash:: build (DistData** inputData, int numItems, int numParents)
  //
  {
    int i = 0;
    int loc = 0;
    int temp = 0;

    // If the data set is empty, then abort.

    if ((numItems <= 0) || (inputData == NULL))
    {
      if (verbosity > 0)
      {
        if (numItems == 1)
        {
          printf ("ERROR (from build): data set has only 1 item.\n");
        }
        else
        {
          printf ("ERROR (from build): empty data set.\n");
        }
        
        fflush (NULL);
      }

      return 0;
    }

    if (verbosity >= 2)
    {
      printf ("Building SASH from data array...\n");
      fflush (NULL);
    }

    data = inputData;

    // Reserve SASH storage, and set up tree parameters.
    // As a result of this operation, the SASH size, number of levels,
    //   etc, are set.

    reserveStorage (numItems, numParents);

    // Randomly assign data items to SASH nodes.

    for (i=0; i<size; i++)
    {
      internToExternMapping[i] = i;
    }

    for (i=size-1; i>=0; i--)
    {
      loc = genInt () % (i+1);
      temp = internToExternMapping[loc];
      internToExternMapping[loc] = internToExternMapping[i];
      internToExternMapping[i] = temp;
    }

    // Recursively build the SASH structure.

    numDistComps = 0UL;

    doBuild (numItems);

    if (verbosity >= 2)
    {
      printStats ();
    }

    return size;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Loads a previously-computed SASH from the specified file.
   * The original data set must also be provided (as well as the
   *   number of items in the data set).
   * The extension ".sash" is automatically appended to the file name.
   * If successful, the number of SASH items is returned.
   * If unsuccessful, zero is returned.
   */

  int Sash:: build (char* fileName, DistData** inputData, int numItems)
  //
  {
    int i = 0;
    int j = 0;
    int loc = 0;
    int temp = 0;
    int inSize = 0;
    int inLevels = 0;
    int inMaxParents = 0;
    int numChildren = 0;
    int* childList = NULL;
    FILE* inFile = NULL;

    // If the data set is empty, then abort.

    if ((fileName == NULL) || (numItems <= 0) || (inputData == NULL))
    {
      if (verbosity > 0)
      {
        if (numItems == 1)
        {
          printf ("ERROR (from build): data set has only 1 item.\n");
        }
        else
        {
          printf ("ERROR (from build): empty data set or filename.\n");
        }
        
        fflush (NULL);
      }

      return 0;
    }

    if (verbosity >= 2)
    {
      printf ("Loading SASH from file %s.sash ...\n", fileName);
      fflush (NULL);
    }

    data = inputData;

    // Open the file containing the SASH.
    // If we fail to open the file, then abort.

    strcpy (stringBuf, fileName);
    strcat (stringBuf, ".sash");
    inFile = fopen (stringBuf, "rt");

    if (inFile == NULL)
    {
      if (verbosity > 0)
      {
        printf
          ("ERROR (from build): file %s.sash could not be opened.\n",
            fileName);
        fflush (NULL);
      }

      return 0;
    }

    // Skip two comment lines.

    fgets (stringBuf, SASH_BUFSIZE_, inFile);
    fgets (stringBuf, SASH_BUFSIZE_, inFile);

    // Read in basic parameters.

    fscanf
        (inFile, "%d%d%d%d%ld",
         &(inSize), &(inLevels), &(inMaxParents), &(numOrphans), &(seed));

    // Are these parameter values what we expected?
    // If not, then abort!

    if (inSize != numItems)
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from build):");
        printf (" unexpected SASH parameters in file %s.sash.\n", fileName);
      }

      fclose (inFile);
      fflush (NULL);
      
      return 0;
    }

    // Reserve SASH storage, and set up tree parameters.
    // As a result of this operation, the expected SASH size,
    //   number of levels, etc, are set.

    reserveStorage (inSize, inMaxParents);

    // Skip another comment line.

    fgets (stringBuf, SASH_BUFSIZE_, inFile);
    fgets (stringBuf, SASH_BUFSIZE_, inFile);

    // Read in information for each node:
    //   internal index,
    //   external index,
    //   number of children,
    //   and indices of children.
    // Build the list of children, if any exist.

    for (i=0; i<size; i++)
    {
      fscanf (inFile, "%d", &(loc));

      if (loc != i)
      {
        if (verbosity > 0)
        {
          printf ("ERROR (from build):");
          printf (" invalid entry in file %s.sash.\n", fileName);
        }
        
        fclose (inFile);
        fflush (NULL);

        return 0;
      }

      fscanf (inFile, "%d%d",
              &(internToExternMapping[i]),
              &(numChildren));

      if (numChildren > 0)
      {
        childList = new int [numChildren];
      }
      else
      {
        childList = NULL;
      }

      childIndexLList[i] = childList;
      childLSizeList[i] = numChildren;

      for (j=0; j<numChildren; j++)
      {
        fscanf (inFile, "%d", &(childList[j]));
      }
    }

    fclose (inFile);
    fflush (NULL);

    if (verbosity >= 2)
    {
      printStats ();
    }

    return size;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Perform an approximate range query for the specified item.
   * The upper limit on the query-to-item distance must be supplied.
   * The number of elements actually found is returned.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findAllInRange (DistData* query, float limit)
  // 
  {
    return findAllInRange (query, limit, 0);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Perform an approximate range query for the specified item.
   * The upper limit on the query-to-item distance must be supplied.
   * The number of elements actually found is returned.
   * The search is relative to a data sample of size N / 2^r,
   *   where N is the number of items in the set, and r is 
   *   a non-negative integer ("sampleRate").
   * A "sampleRate" of zero indicates a search relative to the entire set.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findAllInRange (DistData* query, float limit, int sampleRate)
  // 
  {
    queryResultSize = 0;
    queryResultSampleSize = 0;
    numDistComps = 0UL;

    if ((size <= 0)
        || (query == NULL)
        || (limit < 0.0F)
        || (sampleRate < 0)
        || ((sampleRate >= levels) && (size > 1)))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from findAllInRange): invalid argument(s).\n");
        fflush (NULL);
      }
        
      return 0;
    }

    setNewQuery (query);

    return doFindAllInRange (limit, sampleRate);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Perform an approximate range query for the specified item.
   * The upper limit on the query-to-item distance must be supplied.
   * The number of elements actually found is returned.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findMostInRange (DistData* query, float limit)
  // 
  {
    return findMostInRange (query, limit, 0, 1.0F);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Perform an approximate range query for the specified item.
   * The upper limit on the query-to-item distance must be supplied.
   * The number of elements actually found is returned.
   * The search is relative to a data sample of size N / 2^r,
   *   where N is the number of items in the set, and r is 
   *   a non-negative integer ("sampleRate").
   * A "sampleRate" of zero indicates a search relative to the entire set.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findMostInRange (DistData* query, float limit, int sampleRate)
  // 
  {
    return findMostInRange (query, limit, sampleRate, 1.0F);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Perform an approximate range query for the specified item.
   * The upper limit on the query-to-item distance must be supplied.
   * The number of elements actually found is returned.
   * The method also makes use of a parameter ("scaleFactor")
   *   that influences the trade-off between time and accuracy.
   * The default value of this parameter is 1.0 - increasing the value
   *   will increase running time (roughly proportionally) and increase
   *   the accuracy of the result.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findMostInRange (DistData* query, float limit, float scaleFactor)
  // 
  {
    return findMostInRange (query, limit, 0, scaleFactor);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Perform an approximate range query for the specified item.
   * The upper limit on the query-to-item distance must be supplied.
   * The number of elements actually found is returned.
   * The search is relative to a data sample of size N / 2^r,
   *   where N is the number of items in the set, and r is 
   *   a non-negative integer ("sampleRate").
   * A "sampleRate" of zero indicates a search relative to the entire set.
   * The method also makes use of a parameter ("scaleFactor")
   *   that influences the trade-off between time and accuracy.
   * The default value of this parameter is 1.0 - increasing the value
   *   will increase running time (roughly proportionally) and increase
   *   the accuracy of the result.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findMostInRange
             (DistData* query, float limit, int sampleRate, float scaleFactor)
  // 
  {
    queryResultSize = 0;
    queryResultSampleSize = 0;
    numDistComps = 0UL;

    if ((size <= 0)
        || (query == NULL)
        || (limit < 0.0F)
        || (sampleRate < 0)
        || ((sampleRate >= levels) && (size > 1))
        || (scaleFactor <= 0.0F))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from findMostInRange): invalid argument(s).\n");
        fflush (NULL);
      }
        
      return 0;
    }

    setNewQuery (query);

    return doFindMostInRange (limit, sampleRate, scaleFactor);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Find a set of approximate nearest neighbours for the specified
   *   query item.
   * The desired number of elements must be given ("howMany").
   * The number of elements actually found is returned.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findNear (DistData* query, int howMany)
  // 
  {
    return findNear (query, howMany, 0, 1.0F);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Find a set of approximate nearest neighbours for the specified
   *   query item.
   * The search is relative to a data sample of size N / 2^r,
   *   where N is the number of items in the set, and r is 
   *   a non-negative integer ("sampleRate").
   * A "sampleRate" of zero indicates a search relative to the entire set.
   * The desired number of elements must be given ("howMany").
   * The number of elements actually found is returned.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findNear (DistData* query, int howMany, int sampleRate)
  // 
  {
    return findNear (query, howMany, sampleRate, 1.0F);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Find a set of approximate nearest neighbours for the specified
   *   query item.
   * The desired number of elements must be given ("howMany").
   * The number of elements actually found is returned.
   * The method also makes use of a parameter ("scaleFactor")
   *   that influences the trade-off between time and accuracy.
   * The default value of this parameter is 1.0 - increasing the value
   *   will increase running time (roughly proportionally) and increase
   *   the accuracy of the result.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findNear (DistData* query, int howMany, float scaleFactor)
  // 
  {
    return findNear (query, howMany, 0, scaleFactor);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Find a set of approximate nearest neighbours for the specified
   *   query item.
   * The search is relative to a data sample of size N / 2^r,
   *   where N is the number of items in the set, and r is 
   *   a non-negative integer ("sampleRate").
   * A "sampleRate" of zero indicates a search relative to the entire set.
   * The desired number of elements must be given ("howMany").
   * The number of elements actually found is returned.
   * The method also makes use of a parameter ("scaleFactor")
   *   that influences the trade-off between time and accuracy.
   * The default value of this parameter is 1.0 - increasing the value
   *   will increase running time (roughly proportionally) and increase
   *   the accuracy of the result.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findNear
          (DistData* query, int howMany, int sampleRate, float scaleFactor)
  // 
  {
    queryResultSize = 0;
    queryResultSampleSize = 0;
    numDistComps = 0UL;

    if ((size <= 0)
        || (query == NULL)
        || (howMany <= 0)
        || (sampleRate < 0)
        || ((sampleRate >= levels) && (size > 1))
        || (scaleFactor <= 0.0F))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from findNear): invalid argument(s).\n");
        fflush (NULL);
      }
        
      return 0;
    }

    setNewQuery (query);

    return doFindNear (howMany, sampleRate, scaleFactor);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Find a set of exact nearest neighbours for the specified
   *   query item.
   * The desired number of elements must be given ("howMany").
   * The number of elements actually found is returned.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findNearest (DistData* query, int howMany)
  // 
  {
    return findNearest (query, howMany, 0);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Find a set of exact nearest neighbours for the specified
   *   query item.
   * The search is relative to a data sample of size N / 2^r,
   *   where N is the number of items in the set, and r is 
   *   a non-negative integer ("sampleRate").
   * A "sampleRate" of zero indicates a search relative to the entire set.
   * The desired number of elements must be given ("howMany").
   * The number of elements actually found is returned.
   * The query result can be obtained via calls to the following methods:
   *         getResultAcc
   *         getResultDists
   *         getResultDistComps
   *         getResultIndices
   *         getResultNumFound
   *         getResultSampleSize
   * The result items are sorted in increasing order of their distances
   *   to the query.
   */

  int Sash:: findNearest (DistData* query, int howMany, int sampleRate)
  // 
  {
    queryResultSize = 0;
    queryResultSampleSize = 0;
    numDistComps = 0UL;

    if ((size <= 0)
        || (query == NULL)
        || (howMany <= 0)
        || (sampleRate < 0)
        || ((sampleRate >= levels) && (size > 1)))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from findNearest): invalid argument(s).\n");
        fflush (NULL);
      }
        
      return 0;
    }

    setNewQuery (query);

    return doFindNearest (howMany, sampleRate);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns direct access to the SASH input data list.
   */

  DistData** Sash:: getData ()
  {
    return data;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Fills the supplied list with the mapping from external item 
   *   indices to internal SASH indices.
   * If successful, the number of SASH items is returned.
   * If unsuccessful, zero is returned.
   */

  int Sash:: getExternToInternMapping (int* result, int capacity)
  // 
  {
    int i;

    if ((result == NULL) || (capacity < size))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from getExternToInternMapping): ");
        printf  ("result list capacity is too small.\n");
        fflush (NULL);
      }
        
      return 0;
    }

    for (i=0; i<size; i++)
    {
      result[internToExternMapping[i]] = i;
    }

    return size;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Fills the supplied list with the mapping from internal SASH item 
   *   indices to external indices.
   * If successful, the number of SASH items is returned.
   * If unsuccessful, zero is returned.
   */

  int Sash:: getInternToExternMapping (int* result, int capacity)
  // 
  {
    int i;

    if ((result == NULL) || (capacity < size))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from getInternToExternMapping): ");
        printf  ("result list capacity is too small.\n");
        fflush (NULL);
      }
        
      return 0;
    }

    for (i=0; i<size; i++)
    {
      result[i] = internToExternMapping[i];
    }

    return size;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns the upper limit on the number of parents per SASH node.
   */

  int Sash:: getMaxParents ()
  {
    return maxParents;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns the number of data items of the SASH.
   */

  int Sash:: getNumItems ()
  {
    return size;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns the number of sample levels of the SASH.
   */

  int Sash:: getNumLevels ()
  {
    return levels;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns the number of orphan nodes encountered during SASH construction.
   */

  int Sash:: getNumOrphans ()
  //
  {
    return numOrphans;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Computes the recall accuracy of the most recent query result.
   * A list of the exact distances must be provided, sorted
   *   from smallest to largest.
   * The number of exact distances provided determines the size
   *   of the neighbourhood within which the accuracy is assessed.
   * The list must contain at least as many entries as the number of
   *   items found in the query result.
   * If unsuccessful, a negative value is returned.
   */

  float Sash:: getResultAcc (float* exactDistList, int howMany)
  // 
  {
    int i;
    int loc = 0;

    if ((exactDistList == NULL)
        || (howMany < queryResultSize))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from getResultAcc): ");
        printf  ("exact distance list is too small.\n");
        fflush (NULL);
      }
        
      return SASH_UNKNOWN_;
    }

    for (i=0; i<howMany; i++)
    {
      if ((loc < queryResultSize)
          && (queryResultDistList[loc] <= exactDistList[i]))
      {
        loc++;
      }
    }

    return ((float) loc) / howMany;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Fills the supplied list with the query-to-neighbour
   *   distances found in the most recent SASH query.
   * If successful, the number of items found is returned.
   * If unsuccessful, zero is returned.
   */

  int Sash:: getResultDists (float* result, int capacity)
  // 
  {
    int i;

    if ((result == NULL) || (capacity < queryResultSize))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from getResultDists): ");
        printf  ("result list capacity is too small.\n");
        fflush (NULL);
      }
        
      return 0;
    }

    for (i=0; i<queryResultSize; i++)
    {
      result[i] = queryResultDistList[i];
    }

    return queryResultSize;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns the number of distance computations performed during
   *   the most recent SASH operation.
   */

  unsigned long Sash:: getResultDistComps ()
  // 
  {
    return numDistComps;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Fills the supplied list with the (external) indices of the
   *   items found in the most recent SASH query.
   * If successful, the number of items found is returned.
   * If unsuccessful, zero is returned.
   */

  int Sash:: getResultIndices (int* result, int capacity)
  // 
  {
    int i;

    if ((result == NULL) || (capacity < queryResultSize))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from getResultIndices): ");
        printf  ("result list capacity is too small.\n");
        fflush (NULL);
      }
        
      return 0;
    }

    for (i=0; i<queryResultSize; i++)
    {
      result[i] = internToExternMapping[queryResultIndexList[i]];
    }

    return queryResultSize;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns the number of items found in the most recent query.
   */

  int Sash:: getResultNumFound ()
  // 
  {
    return queryResultSize;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns the sample size used in the most recent query.
   */

  int Sash:: getResultSampleSize ()
  // 
  {
    return queryResultSampleSize;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Returns the seed value used for random number generator initialization.
   */

  unsigned long Sash:: getRNGSeed ()
  //
  {
    return seed;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Fills the supplied list with the mapping from external item
   *   indices to internal SASH sample levels.
   * If successful, the number of SASH items is returned.
   * If unsuccessful, zero is returned.
   */

  int Sash:: getSampleAssignment (int* result, int capacity)
  // 
  {
    int i;
    int lvl;

    if ((result == NULL) || (capacity < size))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from getSampleAssignment): ");
        printf  ("result list capacity is too small.\n");
        fflush (NULL);
      }
        
      return 0;
    }

    for (lvl=0; lvl<levels; lvl++)
    {
      for (i=sampleSizeList[lvl+1]; i<sampleSizeList[lvl]; i++)
      {
        result[internToExternMapping[i]] = lvl;
      }
    }

    result[internToExternMapping[0]] = levels;

    return size;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Fills the supplied list with the SASH sample level sizes,
   *   from smallest to largest.
   * The result does not include the "sample" consisting solely of the
   *   SASH root item.
   * If successful, the number of SASH sample levels is returned
   *   (excluding that of the root).
   * If unsuccessful, zero is returned.
   */

  int Sash:: getSampleSizes (int* result, int capacity)
  // 
  {
    int lvl;

    if ((result == NULL) || (capacity < levels))
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from getSampleSizes): ");
        printf  ("result list capacity is too small.\n");
        fflush (NULL);
      }
        
      return 0;
    }

    for (lvl=0; lvl<levels; lvl++)
    {
      result[lvl] = sampleSizeList[lvl];
    }

    return levels;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Resets the current query object to NULL.
   * This has the effect of clearing any saved distances - subsequent
   *   findNear and findNearest operations would be forced to compute
   *   all needed distances from scratch.
   */

  void Sash:: resetQuery ()
  // 
  {
    setNewQuery (NULL);
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Save the SASH to the specified file.
   * The extension ".sash" is automatically appended to the file name.
   * If successful, the number of SASH items is returned.
   * If unsuccessful, zero is returned.
   */

  int Sash:: saveToFile (char* fileName)
  //
  {
    int i;
    int j;
    int numChildren = 0;
    int* childList = NULL;
    FILE* outFile = NULL;

    // If the SASH has not yet been built, abort.

    if (size <= 0)
    {
      return 0;
    }

    // Open the file for writing.
    // If this fails, then abort.

    if (fileName == NULL)
    {
      if (verbosity > 0)
      {
        printf ("ERROR (from saveToFile): output file name is NULL.\n");
        fflush (NULL);
      }

      return 0;
    }

    strcpy (stringBuf, fileName);
    strcat (stringBuf, ".sash");
    outFile = fopen (stringBuf, "wt");

    if (outFile == NULL)
    {
      if (verbosity > 0)
      {
        printf
          ("ERROR (from saveToFile): file %s.sash could not be opened.\n",
            fileName);
        fflush (NULL);
      }

      return 0;
    }

    // Begin writing to the output file.
    // First, write a comment identifying the SASH version and the
    //   output file name.

    fprintf (outFile, "%% SASH %s: %s.sash\n", SASH_VERSION_, fileName);

    // Write the main SASH parameters.

    fprintf (outFile, "%% size levels maxParents numOrphans seed\n");
    fprintf
        (outFile, "%d %d %d %d %ld\n",
         size, levels, maxParents, numOrphans, seed);

    // For each item, write out:
    //   its index,
    //   the index of the item in the original input list,
    //   the number of children of the item,
    //   and a list of the indices of the children.

    fprintf
        (outFile, "%% nodeID origItemID numChildren child0 child1 ...\n");

    for (i=0; i<size; i++)
    {
      numChildren = childLSizeList[i];
      childList = childIndexLList[i];

      fprintf
          (outFile, "%d %d %d", i, internToExternMapping[i], numChildren);

      for (j=0; j<numChildren; j++)
      {
        fprintf (outFile, " %d", childList[j]);
      }
 
      fprintf (outFile, "\n");
    }

    fclose (outFile);
    fflush (NULL);

    return size;
  }


  //////////////////////////////////////////////////////////////////////


  /**
   * Sets the verbosity level for messages.
   * Verbosity of zero or less: no messages produced.
   * Verbosity of 1: error messages only.
   * Verbosity of 2: error and progress messages only.
   * Verbosity of 3 or more: error, progress, and debug messages reported.
   */

  void Sash:: setVerbosity (int verbosity)
  //
  {
    if (verbosity <= 0)
    {
      this->verbosity = 0;
    }
    else if (verbosity >= 3)
    {
      this->verbosity = 3;
    }
    else
    {
      this->verbosity = verbosity;
    }
  }


  //////////////////////////////////////////////////////////////////////
  //                         Private Methods                          //
  //////////////////////////////////////////////////////////////////////


  float Sash:: computeDistFromQuery (int itemIndex)
  //
  // Returns the distance from the current query object to the 
  //   specified data object.
  // If the distance has already been computed and stored,
  //   the stored distance is returned.
  // Otherwise, the distance is computed and stored before returning it.
  //
  {
    if (distFromQueryList[itemIndex] == SASH_UNKNOWN_)
    {
      distFromQueryList[itemIndex]
        = query->distanceTo (data[internToExternMapping[itemIndex]]);
      storedDistIndexList[numStoredDists] = itemIndex;
      numStoredDists++;
      numDistComps++;
    }

    return distFromQueryList[itemIndex];
  }


  //////////////////////////////////////////////////////////////////////


  void Sash:: doBuild (int numItems)
  // 
  // Recursively builds a SASH on items in the first locations of the
  //   scrambled data array.
  // The number of items to be incorporated into the SASH must be specified.
  // 
  {
    int i = 0;
    int j = 0;
    int parent = 0;
    int child = 0;
    int halfSize = 0;
    int quarterSize = 0;
    unsigned char notFound = TRUE;
    int range = 0;
    float* tempDistList = NULL;
    int* tempIndexList = NULL;

    if (numItems <= maxChildren + 1)
    {
      // We have only a small number of items.
      // Build the SASH explicitly, without recursion.
      // Treat the first array item as the root.

      parentLSizeList[0] = 0;
      parentIndexLList[0] = NULL;

      // Explicitly connect all other items as children of the root,
      //   if they exist.
      // Don't bother sorting the children according to distance from
      //   the root.

      if (numItems == 1)
      {
        levels = 0;
        childLSizeList[0] = 0;
        childIndexLList[0] = NULL;
      }
      else
      {
        // Establish connections from root to children.

        levels = 1;
        childLSizeList[0] = numItems - 1;
        childIndexLList[0] = new int [numItems-1];

        for (i=1; i<numItems; i++)
        {
          childIndexLList[0][i-1] = i;
        }
      }

      if (verbosity >= 2)
      {
        printf ("Number of SASH levels constructed: %d\n", levels);
        fflush (NULL);
      }

      return;
    }

    // The number of items is large.
    // Build the SASH recursively on half the items.

    halfSize = (numItems + 1) / 2;
    doBuild (halfSize);

    // We now want to connect the bottom level of the
    //   recursively-constructed SASH to the remainder of the items.

    // For each item in the remainder, generate a set of tentative
    //   parents from the bottom level of the current SASH.
    // Also, temporarily store (in "childLSizeList") the number of times
    //   each node is requested as a parent.

    for (child=halfSize; child<numItems; child++)
    {
      if ((child % 5000 == 4999) && (verbosity >= 2))
      {
        printf ("Inserting item %d (out of %d)...\n", child+1, size);
        fflush (NULL);
      }

      setNewQuery (data[internToExternMapping[child]]);
      doFindParents (maxParents);

      parentLSizeList[child] = queryResultSize;
      parentIndexLList[child] = new int [maxParents];
      parentDistLList[child] = new float [maxParents];

      for (i=0; i<queryResultSize; i++)
      {
        parentIndexLList[child][i] = queryResultIndexList[i];
        parentDistLList[child][i] = queryResultDistList[i];
        childLSizeList[queryResultIndexList[i]]++;
      }
    }

    // For each parent, reserve tentative storage for its child lists.

    if (halfSize <= maxChildren + 1)
    {
      quarterSize = 1;
    }
    else
    {
      quarterSize = (halfSize + 1) / 2;
    }

    for (parent=quarterSize; parent<halfSize; parent++)
    {
      if (childLSizeList[parent] > 0)
      {
        childIndexLList[parent] = new int [childLSizeList[parent]];
        childDistLList[parent] = new float [childLSizeList[parent]];
        childLSizeList[parent] = 0;
      }
    }

    // Construct tentative child lists for each of the parents,
    //   by reversing the child-to-parent edges.
    // Since the child-to-parent edges are no longer needed,
    //   delete them.

    for (child=halfSize; child<numItems; child++)
    {
      for (i=parentLSizeList[child]-1; i>=0; i--)
      {
        parent = parentIndexLList[child][i];
        j = childLSizeList[parent];
        childIndexLList[parent][j] = child;
        childDistLList[parent][j] = parentDistLList[child][i];
        childLSizeList[parent]++;
      }

      if (parentLSizeList[child] > 0)
      {
        delete [] parentIndexLList[child];
        parentIndexLList[child] = NULL;
        delete [] parentDistLList[child];
        parentDistLList[child] = NULL;
        parentLSizeList[child] = 0;
      }
    }

    // If a child list exceeds the length quota, then trim it.
    // The smallest edges are selected, and the distances are cleared.
    // Simultaneously, inform all children of the number of times
    //   their requests have been granted.
    // This number will be temporarily stored in "parentLSizeList".
    // Also, delete the edge distances as we go, since after the
    //   trimming they are no longer needed.

    for (parent=quarterSize; parent<halfSize; parent++)
    {
      if (childLSizeList[parent] > maxChildren)
      {
        tempDistList = childDistLList[parent];
        tempIndexList = childIndexLList[parent];
        childDistLList[parent] = NULL;
        childIndexLList[parent] = new int [maxChildren];

        childLSizeList[parent]
            = partialQuickSort
                (maxChildren,
                 tempDistList, tempIndexList,
                 0, childLSizeList[parent]-1);

        // Connect the parent to its quota of children.
        // Inform the children that another request has been granted.

        for (i=childLSizeList[parent]-1; i>=0; i--)
        {
          child = tempIndexList[i];
          childIndexLList[parent][i] = child;
          parentLSizeList[child]++;
        }

        delete [] tempDistList;
        tempDistList = NULL;
        delete [] tempIndexList;
        tempIndexList = NULL;
      }
      else
      {
        // Inform the children that another request has been granted.

        for (i=childLSizeList[parent]-1; i>=0; i--)
        {
          parentLSizeList[childIndexLList[parent][i]]++;
        }
      }

      // Eliminate the edge distance information.

      if (childDistLList[parent] != NULL)
      {
        delete [] childDistLList[parent];
        childDistLList[parent] = NULL;
      }
    }

    // For each child, check to see if at least one parent granted
    //   its connection request.
    // (The number of connections granted is stored in "parentLSizeList".)
    // Any "orphans" discovered must be connected to a "foster parent".

    for (child=halfSize; child<numItems; child++)
    {
      if (parentLSizeList[child] == 0)
      {
        // The current child is an orphan.
        // Look for eligible parents by successively widening the range.
        // Eventually a foster parent must be found.
        // But just to be sure, we test to make sure that the range is
        //   not bigger than the number of items in the SASH.

        numOrphans++;
        notFound = TRUE;
        range = 2 * maxParents;

        while (notFound && (range <= numItems))
        {
          setNewQuery (data[internToExternMapping[child]]);
          doFindParents (range);

          for (i=0; i<queryResultSize; i++)
          {
            // Fetch a new candidate foster parent from the query result.

            parent = queryResultIndexList[i];

            // Does this parent have room for another child?
            // If so, then accept the child immediately.

            if (childLSizeList[parent] < maxChildren)
            {
              // Since "childDistLList" is not being used to hold
              //   edge distances any more, we can reuse it to 
              //   indicate whether parents are fostering any orphans.
              // If "childDistLList[parent]!=NULL", then "parent" is
              //   assumed to be fostering at least one orphan.

              if (childDistLList[parent] == NULL)
              {
                // This parent is fostering an orphan for the first time.
                // Expand the size of its child list to the maximum
                //   possible, to accommodate the current orphan
                //   and any future orphans.

                tempIndexList = childIndexLList[parent];
                childIndexLList[parent] = new int [maxChildren];

                for (j=childLSizeList[parent]-1; j>=0; j--)
                {
                  childIndexLList[parent][j] = tempIndexList[j];
                }

                delete [] tempIndexList;
                tempIndexList = NULL;
              }

              // Add the child to the parent's list.
              // To indicate that the parent is now fostering orphans,
              //   set its child edge distance list to anything non-null
              //   (the list "distFromQueryList").

              childDistLList[parent] = distFromQueryList;
              childIndexLList[parent][childLSizeList[parent]] = child;
              childLSizeList[parent]++;
              notFound = FALSE;

              break;
            }
          }

          range *= 2;
        }
      }
    }
    
    // All orphans have now found foster parents.
    // The SASH has grown by one level.
    // Clean up and return.

    for (parent=quarterSize; parent<halfSize; parent++)
    {
      childDistLList[parent] = NULL;
    }

    levels++;

    if (verbosity >= 2)
    {
      printf ("Number of SASH levels constructed: %d\n", levels);
      fflush (NULL);
    }
  }


#ifdef INCREMENTAL_SASH_V1
  int Sash:: addNewChildrens(DistData** newInputData, int newNumItems) {
    int i = 0;
    int loc = 0;
    int temp = 0;

   // TODO: first check if we have enough room to store the new datas

   

    // If the data set is empty, then abort.
    if ((newNumItems <= 0) || (newInputData == NULL)) {
      if (verbosity > 0)
      {
        if (newNumItems == 1)
        {
          printf ("WARNING (from addNewChildrens): data set has only 1 item.\n");
        }
        else
        {
          printf ("WARNING (from addNewChildrens): empty data set.\n");
        }
        
        fflush (NULL);
      }

      return 0;
    }

    if (verbosity >= 2)
    {
      printf ("Adding new datas to SASH from data array...\n");
      fflush (NULL);
    }


    // TODO : HERE we should expand the data
    data = newInputData;

    // Reserve SASH storage, and set up tree parameters.
    // As a result of this operation, the SASH size, number of levels,
    //   etc, are set.

    // TODO: here we should expand the storage 
    reserveStorage (numItems, numParents);



    // Randomly assign data items to SASH nodes.

    for (i=oldsize; i<newsize; i++)
    {
      internToExternMapping[i] = i;
    }

    for (i=newsize-1; i>=oldsize; i--)
    {
      loc = genInt () % (i+1);
      temp = internToExternMapping[loc];
      internToExternMapping[loc] = internToExternMapping[i];
      internToExternMapping[i] = temp;
    }
    // Recursively build the SASH structure.

    numDistComps = 0UL;

    //



    for (child=halfSize; child<numItems; child++)
    {
      if ((child % 5000 == 4999) && (verbosity >= 2))
      {
        printf ("Inserting item %d (out of %d)...\n", child+1, size);
        fflush (NULL);
      }

      setNewQuery (data[internToExternMapping[child]]);
      doFindParents (maxParents);

      parentLSizeList[child] = queryResultSize;
      parentIndexLList[child] = new int [maxParents];
      parentDistLList[child] = new float [maxParents];

      for (i=0; i<queryResultSize; i++)
      {
        parentIndexLList[child][i] = queryResultIndexList[i];
        parentDistLList[child][i] = queryResultDistList[i];
        childLSizeList[queryResultIndexList[i]]++;
      }
    }

    // For each parent, reserve tentative storage for its child lists.

    if (halfSize <= maxChildren + 1)
    {
      quarterSize = 1;
    }
    else
    {
      quarterSize = (halfSize + 1) / 2;
    }

    for (parent=quarterSize; parent<halfSize; parent++)
    {
      if (childLSizeList[parent] > 0)
      {
        childIndexLList[parent] = new int [childLSizeList[parent]];
        childDistLList[parent] = new float [childLSizeList[parent]];
        childLSizeList[parent] = 0;
      }
    }

    // Construct tentative child lists for each of the parents,
    //   by reversing the child-to-parent edges.
    // Since the child-to-parent edges are no longer needed,
    //   delete them.

    for (child=halfSize; child<numItems; child++)
    {
      for (i=parentLSizeList[child]-1; i>=0; i--)
      {
        parent = parentIndexLList[child][i];
        j = childLSizeList[parent];
        childIndexLList[parent][j] = child;
        childDistLList[parent][j] = parentDistLList[child][i];
        childLSizeList[parent]++;
      }

      if (parentLSizeList[child] > 0)
      {
        delete [] parentIndexLList[child];
        parentIndexLList[child] = NULL;
        delete [] parentDistLList[child];
        parentDistLList[child] = NULL;
        parentLSizeList[child] = 0;
      }
    }

    // If a child list exceeds the length quota, then trim it.
    // The smallest edges are selected, and the distances are cleared.
    // Simultaneously, inform all children of the number of times
    //   their requests have been granted.
    // This number will be temporarily stored in "parentLSizeList".
    // Also, delete the edge distances as we go, since after the
    //   trimming they are no longer needed.

    for (parent=quarterSize; parent<halfSize; parent++)
    {
      if (childLSizeList[parent] > maxChildren)
      {
        tempDistList = childDistLList[parent];
        tempIndexList = childIndexLList[parent];
        childDistLList[parent] = NULL;
        childIndexLList[parent] = new int [maxChildren];

        childLSizeList[parent]
            = partialQuickSort
                (maxChildren,
                 tempDistList, tempIndexList,
                 0, childLSizeList[parent]-1);

        // Connect the parent to its quota of children.
        // Inform the children that another request has been granted.

        for (i=childLSizeList[parent]-1; i>=0; i--)
        {
          child = tempIndexList[i];
          childIndexLList[parent][i] = child;
          parentLSizeList[child]++;
        }

        delete [] tempDistList;
        tempDistList = NULL;
        delete [] tempIndexList;
        tempIndexList = NULL;
      }
      else
      {
        // Inform the children that another request has been granted.

        for (i=childLSizeList[parent]-1; i>=0; i--)
        {
          parentLSizeList[childIndexLList[parent][i]]++;
        }
      }

      // Eliminate the edge distance information.

      if (childDistLList[parent] != NULL)
      {
        delete [] childDistLList[parent];
        childDistLList[parent] = NULL;
      }
    }

    // For each child, check to see if at least one parent granted
    //   its connection request.
    // (The number of connections granted is stored in "parentLSizeList".)
    // Any "orphans" discovered must be connected to a "foster parent".

    for (child=halfSize; child<numItems; child++)
    {
      if (parentLSizeList[child] == 0)
      {
        // The current child is an orphan.
        // Look for eligible parents by successively widening the range.
        // Eventually a foster parent must be found.
        // But just to be sure, we test to make sure that the range is
        //   not bigger than the number of items in the SASH.

        numOrphans++;
        notFound = TRUE;
        range = 2 * maxParents;

        while (notFound && (range <= numItems))
        {
          setNewQuery (data[internToExternMapping[child]]);
          doFindParents (range);

          for (i=0; i<queryResultSize; i++)
          {
            // Fetch a new candidate foster parent from the query result.

            parent = queryResultIndexList[i];

            // Does this parent have room for another child?
            // If so, then accept the child immediately.

            if (childLSizeList[parent] < maxChildren)
            {
              // Since "childDistLList" is not being used to hold
              //   edge distances any more, we can reuse it to 
              //   indicate whether parents are fostering any orphans.
              // If "childDistLList[parent]!=NULL", then "parent" is
              //   assumed to be fostering at least one orphan.

              if (childDistLList[parent] == NULL)
              {
                // This parent is fostering an orphan for the first time.
                // Expand the size of its child list to the maximum
                //   possible, to accommodate the current orphan
                //   and any future orphans.

                tempIndexList = childIndexLList[parent];
                childIndexLList[parent] = new int [maxChildren];

                for (j=childLSizeList[parent]-1; j>=0; j--)
                {
                  childIndexLList[parent][j] = tempIndexList[j];
                }

                delete [] tempIndexList;
                tempIndexList = NULL;
              }

              // Add the child to the parent's list.
              // To indicate that the parent is now fostering orphans,
              //   set its child edge distance list to anything non-null
              //   (the list "distFromQueryList").

              childDistLList[parent] = distFromQueryList;
              childIndexLList[parent][childLSizeList[parent]] = child;
              childLSizeList[parent]++;
              notFound = FALSE;

              break;
            }
          }

          range *= 2;
        }
      }
    }
    
    // All orphans have now found foster parents.
    // The SASH has grown by one level.
    // Clean up and return.

    for (parent=quarterSize; parent<halfSize; parent++)
    {
      childDistLList[parent] = NULL;
    }

  }



#endif



  //////////////////////////////////////////////////////////////////////


  int Sash:: doFindAllInRange (float limit, int sampleRate)
  // 
  // Performs an exact range query from the current query object,
  //   with respect to a subset of the items.
  // The subset consists of all items at the indicated sample level and higher.
  // The upper limit on the query-to-item distance is "limit";
  //   the number of neighbours actually found is returned.
  // The results are stored in the SASH query result lists.
  // 
  {
    int i;

    // Handle the singleton case separately.

    if (size == 1)
    {
      queryResultDistList[0] = computeDistFromQuery (0);
      queryResultIndexList[0] = 0;
      queryResultSampleSize = 1;

      if (queryResultDistList[0] <= limit)
      {
        queryResultSize = 1;
      }
      else
      {
        queryResultSize = 0;
      }

      return queryResultSize;
    }

    queryResultSampleSize = sampleSizeList[sampleRate];

    // Compute distances from the current query to all items.

    for (i=0; i<queryResultSampleSize; i++)
    {
      queryResultDistList[i] = computeDistFromQuery (i);
      queryResultIndexList[i] = i;
    }

    // Sort the items by distances, returning the number of 
    //   elements actually found.

    queryResultSize = partialQuickSort
                          (queryResultSampleSize, 
                           queryResultDistList, queryResultIndexList,
                           0, queryResultSampleSize-1);

    // Report only those items whose distances fall within the limit.

    i = 0;
    
    while ((i < queryResultSize) && (queryResultDistList[i] <= limit))
    {
      i++;
    }

    queryResultSize = i;

    return queryResultSize;
  }


  //////////////////////////////////////////////////////////////////////


  int Sash:: doFindMostInRange (float limit, int sampleRate, float scaleFactor)
  // 
  // Performs an approximate range query from the current query object,
  //   with respect to a subset of the items.
  // The subset consists of all items at the indicated sample level and higher.
  // The upper limit on the query-to-item distance must be supplied.
  // The number of elements actually found is returned.
  // The results are stored in the SASH query result lists.
  // The parameter "scaleFactor" influences the tradeoff between speed
  //   and accuracy.
  // The base setting is "scaleFactor=1.0"; queries with "scaleFactor=2.0"
  //   would be expected to be more accurate than those with "scaleFactor=1.0",
  //   but would take roughly twice as much time to process.
  // 
  {
    int i;
    int j;
    int loc;
    int lvl;
    int minNeighbours = 0;
    int activeLevelFirst = 0;
    int activeLevelLast = 0;
    int activeLevelNext = 0;
    int nodeIndex = 0;
    int numChildren = 0;
    int numFound = 0;
    int numSought = 0;

    // Handle the singleton case separately.

    if (size == 1)
    {
      queryResultDistList[0] = computeDistFromQuery (0);
      queryResultIndexList[0] = 0;
      queryResultSampleSize = 1;

      if (queryResultDistList[0] <= limit)
      {
        queryResultSize = 1;
      }
      else
      {
        queryResultSize = 0;
      }

      return queryResultSize;
    }

    // Compute the sample size for the operation.

    queryResultSampleSize = sampleSizeList[sampleRate];

    // Compute the minimum number of neighbours for each sample level.

    minNeighbours
        = (int) ((maxParents * maxChildren * 0.5F * scaleFactor) + 0.5F);

    // Load the root as the tentative sole member of the query result list.
    // If its distance to the query is less than or equal to the limit,
    //   then ensure that it is not overwritten by items from the next
    //   sample level.

    queryResultDistList[0] = computeDistFromQuery (0);
    queryResultIndexList[0] = 0;
    queryResultSize = 1;

    if (queryResultDistList[0] <= limit)
    {
      activeLevelNext = 1;
    }
    else
    {
      activeLevelNext = 0;
    }

    // From the root, search out other nodes to place in the query result.

    for (lvl=levels-1; lvl>=sampleRate; lvl--)
    {
      // For every node at the active level, load its children
      //   into the scratch list, and compute their distances to the query.
      // Nodes in the range [activeLevelFirst..activeLevelNext-1]
      //   have query distances within the limit, and
      //   nodes in the range [activeLevelNext..activeLevelLast]
      //   have query distances in excess of the limit.
      // Either or both of these ranges can be empty.

      scratchListSize = 0;

      for (i=activeLevelFirst; i<=activeLevelLast; i++)
      {
        nodeIndex = queryResultIndexList[i];
        numChildren = childLSizeList[nodeIndex];

        for (j=0; j<numChildren; j++)
        {
          scratchIndexList[scratchListSize] = childIndexLList[nodeIndex][j];
          scratchDistList[scratchListSize]
              = computeDistFromQuery (scratchIndexList[scratchListSize]);
          scratchListSize++;
        }
      }

      // Sort the source lists in place, according to distance.
      // The requested number of edges with smallest distances are preserved,
      //   but other entries may be destroyed.

      numFound = partialQuickSort
                 (scratchListSize, scratchDistList, scratchIndexList,
                  0, scratchListSize-1);

      // Copy over the extracted edges to the output lists,
      //   and return the number of edges extracted.
    
      activeLevelFirst = activeLevelNext;
      loc = 0;

      while ((loc < numFound) && (scratchDistList[loc] <= limit))
      {
        queryResultDistList[activeLevelNext] = scratchDistList[loc];
        queryResultIndexList[activeLevelNext] = scratchIndexList[loc];
        activeLevelNext++;
        loc++;
      }

      // Adjust the set of items for which children will be sought
      //   at the next level.
      // The set will be expanded according to the value of scaleFactor,
      //   provided that this size is at least a certain minimum, and
      //   at most the number of items available.

      activeLevelLast = activeLevelNext - 1;

      numSought = (int) ((loc * scaleFactor) + 0.5F);

      if (numSought < minNeighbours)
      {
        numSought = minNeighbours;
      }

      if (numSought > numFound)
      {
        numSought = numFound;
      }

      while (loc < numSought)
      {
        activeLevelLast++;
        queryResultDistList[activeLevelLast] = scratchDistList[loc];
        queryResultIndexList[activeLevelLast] = scratchIndexList[loc];
        loc++;
      }
    }

    // Sort those items within the range by their distances,
    //   returning the number of elements actually found.

    queryResultSize = partialQuickSort
                          (activeLevelNext,
                           queryResultDistList, queryResultIndexList,
                           0, activeLevelNext-1);

    return queryResultSize;
  }


  //////////////////////////////////////////////////////////////////////


  int Sash:: doFindNear (int howMany, int sampleRate, float scaleFactor)
  // 
  // Computes approximate nearest neighbours of the current query object,
  //   with respect to a subset of the items.
  // The subset consists of all items at the indicated sample level and higher.
  // The number of neighbours sought is "howMany"; the number of neighbours
  //   actually found is returned.
  // The results are stored in the SASH query result lists.
  // The parameter "scaleFactor" influences the tradeoff between speed
  //   and accuracy.
  // The base setting is "scaleFactor=1.0"; queries with "scaleFactor=2.0"
  //   would be expected to be more accurate than those with "scaleFactor=1.0",
  //   but would take roughly twice as much time to process.
  // 
  {
    int i;
    int j;
    int lvl;
    int minNeighbours = 0;
    int levelQuota = 0;
    double reductionFactor = 0.0F;
    int activeLevelFirst = 0;
    int activeLevelLast = 0;
    int nodeIndex = 0;
    int numChildren = 0;
    int numFound = 0;

    // Handle the singleton case separately.

    if (size == 1)
    {
      queryResultSize = 1;
      queryResultDistList[0] = computeDistFromQuery (0);
      queryResultIndexList[0] = 0;
      queryResultSampleSize = 1;

      return 1;
    }

    // Compute the sample size for the operation.

    queryResultSampleSize = sampleSizeList[sampleRate];

    // Compute the item quota for each sample level.

    minNeighbours
        = (int) ((maxParents * maxChildren * 0.5F * scaleFactor) + 0.5F);

    reductionFactor = 1.0F
                      /
                      pow (
                            (double) howMany,
                            1.0F
                            /
                            (
                              ( 
                                log ((double) size)
                                /
                                log (2.0F)
                              )
                              -
                              sampleRate
                            )
                          );

    levelQuota = (int) ((howMany * scaleFactor) + 0.5F);

    for (lvl=sampleRate; lvl<levels; lvl++)
    {
      if (levelQuota < minNeighbours)
      {
        levelQuota = minNeighbours;
      }

      levelQuotaList[lvl] = levelQuota;

      levelQuota = (int) ((reductionFactor * levelQuota) + 0.5F);
    }

    // Load the root as the tentative sole member of the query result list.

    queryResultDistList[0] = computeDistFromQuery (0);
    queryResultIndexList[0] = 0;
    queryResultSize = 1;

    // From the root, search out other nodes to place in the query result.

    for (lvl=levels-1; lvl>=sampleRate; lvl--)
    {
      // For every node at the active level, load its children
      //   into the scratch list, and compute their distances to the query.

      scratchListSize = 0;

      for (i=activeLevelFirst; i<=activeLevelLast; i++)
      {
        nodeIndex = queryResultIndexList[i];
        numChildren = childLSizeList[nodeIndex];

        for (j=0; j<numChildren; j++)
        {
          scratchIndexList[scratchListSize] = childIndexLList[nodeIndex][j];
          scratchDistList[scratchListSize]
              = computeDistFromQuery (scratchIndexList[scratchListSize]);
          scratchListSize++;
        }
      }

      // Extract the closest nodes from the list of accumulated children,
      //   and append them to the tentative query result.

      numFound = extractBestEdges
                      (levelQuotaList[lvl],
                       queryResultDistList, queryResultIndexList,
                       activeLevelLast+1, size,
                       scratchDistList, scratchIndexList,
                       0, scratchListSize-1);

      activeLevelFirst = activeLevelLast + 1;
      activeLevelLast += numFound;
    }

    queryResultSize = activeLevelLast + 1;

    // Sort the items by distances, returning the number of 
    //   elements actually found.

    queryResultSize = partialQuickSort
                          (howMany,
                           queryResultDistList, queryResultIndexList,
                           0, queryResultSize-1);

    return queryResultSize;
  }


  //////////////////////////////////////////////////////////////////////


  int Sash:: doFindNearest (int howMany, int sampleRate)
  // 
  // Computes exact nearest neighbours of the current query object,
  //   with respect to a subset of the items.
  // The subset consists of all items at the indicated sample level and higher.
  // The number of neighbours sought is "howMany"; the number of neighbours
  //   actually found is returned.
  // The results are stored in the SASH query result lists.
  // 
  {
    int i;

    // Handle the singleton case separately.

    if (size == 1)
    {
      queryResultSize = 1;
      queryResultDistList[0] = computeDistFromQuery (0);
      queryResultIndexList[0] = 0;
      queryResultSampleSize = 1;

      return 1;
    }

    queryResultSampleSize = sampleSizeList[sampleRate];

    // Compute distances from the current query to all items.

    for (i=0; i<queryResultSampleSize; i++)
    {
      queryResultDistList[i] = computeDistFromQuery (i);
      queryResultIndexList[i] = i;
    }

    // Sort the items by distances, returning the number of 
    //   elements actually found.

    queryResultSize = partialQuickSort
                          (howMany,
                           queryResultDistList, queryResultIndexList,
                           0, queryResultSampleSize-1);

    return queryResultSize;
  }


  //////////////////////////////////////////////////////////////////////


  int Sash:: doFindParents (int howMany)
  // 
  // Finds a set of parents for the current query item from among the
  //   bottom-level items of the current SASH.
  // The results are stored in the SASH query result lists.
  // 
  {
    int i;
    int j;
    int lvl;
    int nodeIndex = 0;
    int numChildren = 0;

    // Load the root as the tentative sole member of the query result list.

    queryResultDistList[0] = computeDistFromQuery (0);
    queryResultIndexList[0] = 0;
    queryResultSize = 1;

    // From the root, search out other nodes to place in the query result.

    for (lvl=0; lvl<levels; lvl++)
    {
      // For every node at the active level, load its children
      //   into the scratch list, and compute their distances to the query.

      scratchListSize = 0;

      for (i=0; i<queryResultSize; i++)
      {
        nodeIndex = queryResultIndexList[i];
        numChildren = childLSizeList[nodeIndex];

        for (j=0; j<numChildren; j++)
        {
          scratchIndexList[scratchListSize] = childIndexLList[nodeIndex][j];
          scratchDistList[scratchListSize]
              = computeDistFromQuery (scratchIndexList[scratchListSize]);
          scratchListSize++;
        }
      }

      // Extract the closest nodes from the list of accumulated children,
      //   and keep them as the tentative parents of the query.
      // Note that only the candidates from the most recently-processed
      //   level are kept.

      queryResultSize = extractBestEdges
                              (howMany,
                               queryResultDistList, queryResultIndexList,
                               0, size,
                               scratchDistList, scratchIndexList,
                               0, scratchListSize-1);
    }

    return queryResultSize;
  }


  //////////////////////////////////////////////////////////////////////


  int Sash:: extractBestEdges
      (int howMany,
       float* toDistList, int* toIndexList, int toFirst, int toCapacity,
       float* fromDistList, int* fromIndexList, int fromFirst, int fromLast)
  // 
  // Copies a requested number of directed edges having minimum distances
  //   to their targets.
  // The input edges are stored in "fromIndexList" and "fromDistList", in the
  //   range of locations beginning at "fromFirst" and ending at "fromLast".
  // The extracted edges are stored in "toIndexList" and "toDistList",
  //   in locations starting at "toFirst".
  // If the requested number of edges "howMany" would exceed the 
  //   output list capacity "toCapacity" if copying were to start at
  //   location "toFirst", then the number of extracted edges is reduced.
  // WARNING: this operation destroys entries of the input list!
  // 
  {
    int i;
    int numExtracted = 0;

    // Make sure that we don't extract more edges than currently exist.
    
    if (howMany > toCapacity - toFirst)
    {
      howMany = toCapacity - toFirst;
    }

    // Sort the source lists in place, according to distance.
    // The requested number of edges with smallest distances are preserved,
    //   but other entries may be destroyed.

    numExtracted = partialQuickSort
                   (howMany, fromDistList, fromIndexList, fromFirst, fromLast);

    // Copy over the extracted edges to the output lists,
    //   and return the number of edges extracted.
    
    for (i=0; i<numExtracted; i++)
    {
      toDistList[toFirst] = fromDistList[fromFirst];
      toIndexList[toFirst] = fromIndexList[fromFirst];
      toFirst++;
      fromFirst++;
    }

    return numExtracted;
  }


  //////////////////////////////////////////////////////////////////////


  int Sash:: partialQuickSort
              (int howMany,
               float* distList, int* indexList,
               int rangeFirst, int rangeLast)
  //
  // Sorts the smallest items in the supplied list ranges, in place,
  //   according to distances.
  // A partial quicksort is used to sort only the requested number
  //   of items.
  // The smallest items are placed at the beginning of the range, 
  //   in increasing order of distance.
  // WARNING: the remainder of the range can become corrupted
  //   by this operation!
  //
  {
    int i;
    int pivotLoc = 0;
    int pivotIndex = 0;
    float pivotDist = 0.0F;
    int tempIndex = 0;
    float tempDist = 0.0F;
    int low = 0;
    int high = 0;
    int numFound = 0;
    int numDuplicatesToReplace = 0;
    int tieBreakIndex = 0;

    // If the range is empty, or if we've been asked to sort no
    //   items, then return immediately.

    if ((rangeLast < rangeFirst) || (howMany < 1))
    {
      return 0;
    }

    // If there is exactly one element, then again there is nothing
    //   that need be done.

    if (rangeLast == rangeFirst)
    {
      return 1;
    }

    // If the range to be sorted is small, just do an insertion sort.

    if (rangeLast - rangeFirst < 7)
    {
      high = rangeFirst + 1;
      tieBreakIndex = indexList[genInt () % (rangeLast-rangeFirst+1)];

      // The outer while loop considers each item in turn (starting
      //   with the second item in the range), for insertion into
      //   the sorted list of items that precedes it.

      while (high <= rangeLast)
      {
        // Copy the next item to be inserted, as the "pivot".
        // Start the insertion tests with its immediate predecessor.

        pivotDist = distList[high];
        pivotIndex = indexList[high];
        low = high - 1;

        // Work our way down through previously-sorted items
        //   towards the start of the range.

        while (low >= rangeFirst)
        {
          // Compare the item to be inserted (the "pivot") with
          //   the current item.

          if (distList[low] < pivotDist)
          {
            // The current item precedes the pivot in the sorted order.
            // Break out of the loop - we have found the insertion point.

            break;
          }
          else if (distList[low] > pivotDist)
          {
            // The current item follows the pivot in the sorted order.
            // Shift the current item one spot upwards, to make room
            //   for inserting the pivot below it.

            distList[low+1] = distList[low];
            indexList[low+1] = indexList[low];
            low--;
          }
          else
          {
            if (indexList[low] != pivotIndex)
            {
              // The items have the same sort value but are not identical.
              // Break the tie pseudo-randomly.

              if (
                   (
                     (tieBreakIndex < pivotIndex)
                     &&
                     (tieBreakIndex < indexList[low])
                     &&
                     (indexList[low] < pivotIndex)
                   )
                   ||
                   (
                     (tieBreakIndex >= pivotIndex)
                     &&
                     (
                       (tieBreakIndex < indexList[low])
                       ||
                       (indexList[low] < pivotIndex)
                     )
                   )
                 )
              {
                // The current item precedes the pivot in the sorted order.
                // Break out of the loop - we have found the insertion point.

                break;
              }
              else
              {
                // The current item follows the pivot in the sorted order.
                // Shift the current item one spot upwards, to make room
                //   for inserting the pivot below it.

                distList[low+1] = distList[low];
                indexList[low+1] = indexList[low];
                low--;
              }
            }
            else
            {
              // Oh no!
              // We opened up an empty slot for the pivot,
              //   only to find that it's a duplicate of the current item!
              // Close the slot up again, and eliminate the duplicate.

              for (i=low+1; i<high; i++)
              {
                distList[i] = distList[i+1];
                indexList[i] = indexList[i+1];
              }

              // To eliminate the duplicate, overwrite its location with the
              //   item from the end of the range, and then shrink the range
              //   by one.

              distList[high] = distList[rangeLast];
              indexList[high] = indexList[rangeLast];
              rangeLast--;

              // The next iteration must not advance "high", since we've
              //   just put a new element into it which needs to be processed. 
              // Decrementing it here will cancel out with the incrementation
              //   of the next iteration.

              high--;

              // When we break the loop, the pivot element will be put
              //   in its proper place ("low" + 1)
              // Here, the proper place is where rangeLast used to be.
              // To achieve this, we need to adjust "low" here.

              low = rangeLast;

              break;
            }
          }
        }

        // If we've made it to here, we've found the insertion
        //   spot for the current element.
        // Perform the insertion.

        low++;
        distList[low] = pivotDist;
        indexList[low] = pivotIndex;

        // Move to the next item to be inserted in the growing sorted list.

        high++;
      }

      // Return the number of sorted items found.

      numFound = rangeLast - rangeFirst + 1;

      if (numFound > howMany)
      {
        numFound = howMany;
      }

      return numFound;
    }

    // The range to be sorted is large, so do a partial quicksort.
    // Select a pivot item, and swap it with the item at the beginning
    //   of the range.

    pivotLoc = rangeFirst + (genInt () % (rangeLast-rangeFirst+1));
    tieBreakIndex = indexList[genInt () % (rangeLast-rangeFirst+1)];

    pivotDist = distList[pivotLoc];
    distList[pivotLoc] = distList[rangeFirst];
    distList[rangeFirst] = pivotDist;

    pivotIndex = indexList[pivotLoc];
    indexList[pivotLoc] = indexList[rangeFirst];
    indexList[rangeFirst] = pivotIndex;

    // Eliminate all duplicates of the pivot.
    // Any duplicates found are pushed to the end of the range, and
    //   the range shrunk by one (thereby excluding them).

    i = rangeFirst + 1;

    while (i <= rangeLast)
    {
      if ((pivotDist == distList[i]) && (pivotIndex == indexList[i]))
      {
        distList[i] = distList[rangeLast];
        indexList[i] = indexList[rangeLast];
        rangeLast--;
      }
      else
      {
        i++;
      }
    }

    // Partition the remaining items with respect to the pivot.
    // This efficient method is adapted from the one outlined in
    //   Cormen, Leiserson & Rivest.
    // The range is scanned from both ends.
    // Items with small distances are placed below "low", and those
    //   with large distances are placed above "high".
    // Where "low" and "high" meet, the pivot item is inserted.

    low = rangeFirst;
    high = rangeLast + 1;

    while (TRUE)
    {
      // Move the "high" endpoint down until it meets either the pivot,
      //   or something that belongs on the "low" side.
      // If the key values are tied, decide pseudo-randomly.

      do
      {
        high--;
      }
      while (
              (distList[high] > pivotDist)
              ||
              (
                (distList[high] == pivotDist)
                &&
                (
                  (
                    (tieBreakIndex >= pivotIndex)
                    &&
                    (pivotIndex < indexList[high])
                    &&
                    (indexList[high] <= tieBreakIndex)
                  )
                  ||
                  (
                    (tieBreakIndex < pivotIndex)
                    &&
                    (
                      (pivotIndex < indexList[high])
                      ||
                      (indexList[high] <= tieBreakIndex)
                    )
                  )
                )
              )
            );

      // Move the "low" endpoint up until it meets either the pivot,
      //   or something that belongs on the "high" side.
      // If the key values are tied, decide pseudo-randomly.

      do
      {
        low++;
      }
      while (
              (low < high) 
              &&
              (
                (distList[low] < pivotDist)
                ||
                (
                  (distList[low] == pivotDist)
                  &&
                  (
                    (
                      (tieBreakIndex < pivotIndex)
                      &&
                      (tieBreakIndex < indexList[low])
                      &&
                      (indexList[low] < pivotIndex)
                    )
                    ||
                    (
                      (tieBreakIndex >= pivotIndex)
                      &&
                      (
                        (tieBreakIndex < indexList[low])
                        ||
                        (indexList[low] < pivotIndex)
                      )
                    )
                  )
                )
              )
            );

      // Have the "low" and "high" endpoints crossed?
      // If not, we still have more work to do.

      if (low < high)
      {
        // Swap the misplaced items, and try again.

        tempDist = distList[low];
        distList[low] = distList[high];
        distList[high] = tempDist;

        tempIndex = indexList[low];
        indexList[low] = indexList[high];
        indexList[high] = tempIndex;
      }
      else
      {
        // We found the cross-over point.

        break;
      }
    }

    // The pivot value ends up at the location referenced by "high".
    // Swap it with the pivot (which resides at the beginning of the range).

    distList[rangeFirst] = distList[high];
    distList[high] = pivotDist;

    indexList[rangeFirst] = indexList[high];
    indexList[high] = pivotIndex;

    pivotLoc = high;

    // The partition is complete.
    // Recursively sort the items with smaller distance.

    numFound = partialQuickSort
               (howMany, distList, indexList, rangeFirst, pivotLoc-1);

    // If we found enough items (including the pivot), then we are done.
    // Make sure the pivot is in its correct position, if it is used.

    if (numFound >= howMany - 1)
    {
      if (numFound == howMany - 1)
      {
        distList[rangeFirst+numFound] = pivotDist;
        indexList[rangeFirst+numFound] = pivotIndex;
      }

      return howMany;
    }

    // We didn't find enough items, even taking the pivot into account.
    // Were any duplicates discovered during this call?

    if (numFound < pivotLoc - rangeFirst)
    {
      // Duplicates were discovered!
      // Figure out the minimum number of duplicates that must be
      //   replaced by items from the end of the range in order to
      //   leave the non-duplicates in contiguous locations.

      numDuplicatesToReplace = pivotLoc - rangeFirst - numFound;
      high = rangeLast;

      if (numDuplicatesToReplace > rangeLast - pivotLoc)
      {
        numDuplicatesToReplace = rangeLast - pivotLoc;
        rangeLast = rangeFirst + numFound + numDuplicatesToReplace;
      }
      else
      {
        rangeLast -= numDuplicatesToReplace;
      }

      // Replace the required number of duplicates by items from
      //   the end of the range.
      // The size of the range will shrink as a result.

      low = rangeFirst + numFound + 1;

      for (i=0; i<numDuplicatesToReplace; i++)
      {
        distList[low] = distList[high];
        indexList[low] = indexList[high];
        low++;
        high--;
      }
    }

    // Put the pivot element in its proper place.

    distList[rangeFirst+numFound] = pivotDist;
    indexList[rangeFirst+numFound] = pivotIndex;

    // Finish up by sorting larger-distance items.
    // Note that the number of sorted items needed has dropped.

    return numFound + 1 + partialQuickSort
                                (howMany - numFound - 1,
                                 distList, indexList,
                                 rangeFirst+numFound+1, rangeLast);
  }


  //////////////////////////////////////////////////////////////////////


  void Sash:: printStats ()
  //
  // Print statistics related to the SASH construction.
  // Should only be called immediately after the construction.
  //
  {
    printf ("\n");
    printf ("SASH build statistics:\n");
    printf ("  size                  == %d\n", size);
    printf ("  levels                == %d\n", levels);
    printf ("  max parents per node  == %d\n", maxParents);
    printf ("  max children per node == %d\n", maxChildren);
    printf ("  orphan nodes          == %d\n", numOrphans);
    printf ("  distance comparisons  == %ld\n", numDistComps);
    printf ("  RNG seed              == %ld\n", seed);
    printf ("\n");
    fflush (NULL);
  }


  //////////////////////////////////////////////////////////////////////

#if WITH_SASH_INCREMENTAL_V1
  void Sash:: reserveStorage (int numItems, int numParents)
  //
  // Reserve storage for the SASH and its data.
  // The number of SASH items and the maximum number of parents per node
  //   must be given.
  //
  {


    int i;
    int sampleSize = 0;

    // MAKE IT SIMPLE JUST INSTANTIATE TWICE THE AMOUNT OF DATE FOR THE MOMENT
    numItems*=2;

    size = numItems;

    if (numParents <= 3)
    {
      maxParents = 3;
    }
    else
    {
      maxParents = numParents;
    }

    maxChildren = 4 * maxParents;

    // Determine the number of SASH levels, and
    //   the level sample sizes.

    levels = 0;

    if (size > 1)
    {
      levels = 1;
      sampleSize = size;

      while (sampleSize > maxChildren + 1)
      {
        levels++;
        sampleSize = (sampleSize + 1) / 2;
      }
    }

    levelQuotaList = new int [levels];
    sampleSizeList = new int [levels+1];
    sampleSize = size;

    for (i=0; i<levels; i++)
    {
      levelQuotaList[i] = 0;
      sampleSizeList[i] = sampleSize;
      sampleSize = (sampleSize + 1) / 2;
    }

    sampleSizeList[levels] = 1;

    // Reserve storage for the mapping between internal and external
    //   data indices.

    internToExternMapping = new int [size];

    for (i=0; i<size; i++)
    {
      internToExternMapping[i] = i;
    }

    // Set up storage for child-to-parent edges and parent-to-child edges.

    parentIndexLList = new int* [size];
    parentDistLList = new float* [size];
    parentLSizeList = new int [size];

    childIndexLList = new int* [size];
    childDistLList = new float* [size];
    childLSizeList = new int [size];

    for (i=0; i<size; i++)
    {
      parentIndexLList[i] = NULL;
      parentDistLList[i] = NULL;
      parentLSizeList[i] = 0;

      childIndexLList[i] = NULL;
      childDistLList[i] = NULL;
      childLSizeList[i] = 0;
    }

    // Set up storage for managing distance computations and
    //   query results.

    distFromQueryList = new float [size];
    storedDistIndexList = new int [size];
    numStoredDists = 0;

    queryResultDistList = new float [size];
    queryResultIndexList = new int [size];
    queryResultSize = 0;
    queryResultSampleSize = 0;

    for (i=0; i<size; i++)
    {
      distFromQueryList[i] = SASH_UNKNOWN_;
      storedDistIndexList[i] = SASH_NONE_;

      queryResultDistList[i] = SASH_UNKNOWN_;
      queryResultIndexList[i] = SASH_NONE_;
    }

    // Set up temporary lists for edge sorting and accumulation
    //   during searches.

    scratchListSize = maxChildren * ((size + 1) / 2);
    scratchDistList = new float [scratchListSize];
    scratchIndexList = new int [scratchListSize];

    for (i=0; i<scratchListSize; i++)
    {
      scratchDistList[i] = SASH_UNKNOWN_;
      scratchIndexList[i] = SASH_NONE_;
    }

    scratchListSize = 0;
  }
#else
  void Sash:: reserveStorage (int numItems, int numParents)
  //
  // Reserve storage for the SASH and its data.
  // The number of SASH items and the maximum number of parents per node
  //   must be given.
  //
  {
    int i;
    int sampleSize = 0;

    size = numItems;

    if (numParents <= 3)
    {
      maxParents = 3;
    }
    else
    {
      maxParents = numParents;
    }

    maxChildren = 4 * maxParents;

    // Determine the number of SASH levels, and
    //   the level sample sizes.

    levels = 0;

    if (size > 1)
    {
      levels = 1;
      sampleSize = size;

      while (sampleSize > maxChildren + 1)
      {
        levels++;
        sampleSize = (sampleSize + 1) / 2;
      }
    }

    levelQuotaList = new int [levels];
    sampleSizeList = new int [levels+1];
    sampleSize = size;

    for (i=0; i<levels; i++)
    {
      levelQuotaList[i] = 0;
      sampleSizeList[i] = sampleSize;
      sampleSize = (sampleSize + 1) / 2;
    }

    sampleSizeList[levels] = 1;

    // Reserve storage for the mapping between internal and external
    //   data indices.

    internToExternMapping = new int [size];

    for (i=0; i<size; i++)
    {
      internToExternMapping[i] = i;
    }

    // Set up storage for child-to-parent edges and parent-to-child edges.

    parentIndexLList = new int* [size];
    parentDistLList = new float* [size];
    parentLSizeList = new int [size];

    childIndexLList = new int* [size];
    childDistLList = new float* [size];
    childLSizeList = new int [size];

    for (i=0; i<size; i++)
    {
      parentIndexLList[i] = NULL;
      parentDistLList[i] = NULL;
      parentLSizeList[i] = 0;

      childIndexLList[i] = NULL;
      childDistLList[i] = NULL;
      childLSizeList[i] = 0;
    }

    // Set up storage for managing distance computations and
    //   query results.

    distFromQueryList = new float [size];
    storedDistIndexList = new int [size];
    numStoredDists = 0;

    queryResultDistList = new float [size];
    queryResultIndexList = new int [size];
    queryResultSize = 0;
    queryResultSampleSize = 0;

    for (i=0; i<size; i++)
    {
      distFromQueryList[i] = SASH_UNKNOWN_;
      storedDistIndexList[i] = SASH_NONE_;

      queryResultDistList[i] = SASH_UNKNOWN_;
      queryResultIndexList[i] = SASH_NONE_;
    }

    // Set up temporary lists for edge sorting and accumulation
    //   during searches.

    scratchListSize = maxChildren * ((size + 1) / 2);
    scratchDistList = new float [scratchListSize];
    scratchIndexList = new int [scratchListSize];

    for (i=0; i<scratchListSize; i++)
    {
      scratchDistList[i] = SASH_UNKNOWN_;
      scratchIndexList[i] = SASH_NONE_;
    }

    scratchListSize = 0;
  }
#endif



  //////////////////////////////////////////////////////////////////////


  void Sash:: setNewQuery (DistData* query)
  //
  // Accepts a new item as the query object for future distance comparisons.
  // Any previously-stored distances are cleared by this operation,
  //   except in the case where the previous query object is identical
  //   to the current query object.
  //
  {
    int i;

    if (query != this->query)
    {
      for (i=0; i<numStoredDists; i++)
      {
        distFromQueryList[storedDistIndexList[i]] = SASH_UNKNOWN_;
      }

      this->query = query;
      numStoredDists = 0;
    }
  }


  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////

